
;;============================================================
;; Routine to create comparison plots between the various rho_xx being
;; considered in this study.

subset = ''

;; Define the COMMON BLOCK
COMMON FFORE_BLOCK,n,rb3,R0,d,R,Z,tau,f_data,sig_f,corr,farlist,do_Tfit,$
   lpstr,rho_hi,rho_h2,rho_star

FORWARD_FUNCTION BGPS_FFORE_MODEL, BGPS_FFORE_MODEL2, BGPS_FFORE_SIMPLE, $
   STAR_DENSITY, BGPS_LOOPIE

RESOLVE_ROUTINE, 'BGPS_FFORE', /COMPILE_FULL_FILE, /EITHER, /NO_RECOMPILE

restore,'./irdc_dist_model/bgps_rb3'+subset+'.sav',/ver
n = n_elements(rb3)

;; Define arrays for TEMP_DEPENDENT_VARIABLES
tau    = dblarr(n)
f_data = dblarr(n)
sig_f  = dblarr(n)

;; Get galactic params
defsysv, '!MW', exists = exists
IF NOT exists THEN galactic_params 
;; Using !MW parameters -- Goes into COMMON block!
R0 = !MW.R0
d = dindgen(!MW.NBINS*5.)*!MW.BINSIZE + !MW.BINSTART
corr = R0 / 8.5d3

myps,'./irdc_dist_model/analysis_plots/rho_comparisons.eps',xsize=10,ysize=4.5

multiplot,[2,1],xgap=0.04
rhosym = cgSymbol('rho')

;; First panel -- Along plane

R = findgen(1001)*10.
Z = R*0.

plot,R/1.d3,Z,/nodata,yr=[1.d-1,1.d2],/ylog,xtit='R!dgal!n [kpc]  (Z = 0)',$
     ytickformat='exponent10',ytit=rhosym+' / '+rhosym+'(R!do!n,0)'

ind = WHERE( abs(R - R0) EQ min( abs(R - R0) ), nind)

;; Get densities
rho_h2 = h2_density(R,Z)
rho_h2p = rho_h2 / (rho_h2[ind])[0]

rho_hi = hi_density(R,Z)
rho_hip = rho_hi / (rho_hi[ind])[0]

rho_star = star_density(R,Z)
rho_sdp = rho_star / (rho_star[ind])[0]



;; Get fit densities -- so load files

colors = ['PBG3','Orchid','TG3']

FOR model=6,4,-1 DO BEGIN

   restore,'./irdc_dist_model/data/case'+string(model,format="(I0)")+$
           '_irdc_TdFIT.sav'
   
   
   CASE model OF
      4: junk = BGPS_FFORE_MODEL(R0, res[0],res[1],res[2],res[3],rho)
      5: junk = BGPS_FFORE_MODEL2(R0, res[0],res[1],res[2],rho)
      6: junk = BGPS_FFORE_SIMPLE(R0, res[0],res[1],res[2],rho)
   ENDCASE
   
   rho /= (rho[ind])[0]
   
   oplot,R/1.d3,rho,color=cgColor(colors[model-4]),thick=model-1
   
   
ENDFOR

oplot,R/1.d3,rho_h2p,color=cgColor('GRN4')
oplot,R/1.d3,rho_hip,color=cgColor('BLK5')
;oplot,R/1.d3,rho_sdp,color=cgColor('BLK5')


al_legend,/top,/right,box=0,linestyle=0,thick=[2,3,4,1,1],$
          color=[colors,'BLK5','GRN4'],charsize=0.8,$
          ['Case '+string(indgen(3)+4,format="(I0)"),$
           rhosym+'!dHI!n',rhosym+'!dH2!n']

multiplot

;; Second panel -- Vertically

Z = findgen(1001)/2.5 - 200.
R = Z*0. + R0


plot,Z,R,/nodata,yr=[1.d-1,2],/ylog,/xst,/yst,$
     xtit='Z [pc]  (R!dgal!n = R!do!n)',$
     ytickformat='exponent10'

ind = WHERE( Z EQ 0, nind)
print,nind

;; Get densities
rho_h2 = h2_density(R,Z)
rho_h2p = rho_h2 / (rho_h2[ind])[0]

rho_hi = hi_density(R,Z)
rho_hip = rho_hi / (rho_hi[ind])[0]

rho_star = star_density(R,Z)
rho_sdp = rho_star / (rho_star[ind])[0]

;; Get fit densities -- so load files

colors = ['PBG3','Orchid','TG3']

FOR model=6,4,-1 DO BEGIN

   restore,'./irdc_dist_model/data/case'+string(model,format="(I0)")+$
           '_irdc_TdFIT.sav'
   
   
   CASE model OF
      4: junk = BGPS_FFORE_MODEL(R0, res[0],res[1],res[2],res[3],rho)
      5: junk = BGPS_FFORE_MODEL2(R0, res[0],res[1],res[2],rho)
      6: junk = BGPS_FFORE_SIMPLE(R0, res[0],res[1],res[2],rho)
   ENDCASE
   
   rho /= (rho[ind])[0]
   
   oplot,Z,rho,color=cgColor(colors[model-4]),thick=model-1
   
   
ENDFOR

oplot,Z,rho_h2p,color=cgColor('GRN4')
oplot,Z,rho_hip,color=cgColor('BLK5')
;oplot,Z,rho_sdp,color=cgColor('BLK5')


al_legend,/top,/right,box=0,linestyle=0,thick=[2,3,4,1,1],$
          color=[colors,'BLK5','GRN4'],charsize=0.8,$
          ['Case '+string(indgen(3)+4,format="(I0)"),$
           rhosym+'!dHI!n',rhosym+'!dH2!n']


multiplot,/reset
multiplot,/default
multiplot,/reset


myps,/done

END
